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Abstract. 

We present a general methodology for determining the gamma-ray flux from annihilation 
of dark matter particles in Milky Way satellite galaxies, focusing on two promising satellites 
as examples: Segue 1 and Draco. We use the SuperBayeS code to explore the best-fitting 
regions of the Constrained Minimal Supersymmetric Standard Model (CMSSM) parameter 
space, and an independent MCMC analysis of the dark matter halo properties of the satelUtes 
using published radial velocities. We present a formalism for determining the boost from halo 
substructure in these galaxies and show that its value depends strongly on the extrapolation 
of the concentration-mass (c(M)) relation for CDM subhalos down to the minimum possible 
mass. We show that the preferred region for this minimum halo mass within the CMSSM with 
neutralino dark matter is ~ 10"' - 10"^ M©. For the boost model where the observed power- 
law c(M) relation is extrapolated down to the minimum halo mass we find average boosts 
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of about 20, while the Bullock et al (2001) c(M) model results in boosts of order unity. We 
estimate that for the power-law c(M) boost model and photon energies greater than a GeV, the 
Fermi space-telescope has about 20% chance of detecting a dark matter annihilation signal 
from Draco with signal-to-noise greater than 3 after about 5 years of observation. 
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1. Introduction 

If the dark matter in the Universe is comprised of stable, weakly-interacting massive particles 
(WIMPs), in many instances this leads to the prediction that WIMPs will self-annihilate 
into Standard Model (SM) particles that may be visible with the upcoming generation of 
high-energy particle detectors. If high energy gamma-rays are produced, there are several 
promising sources within our own Galactic environment where the annihilation radiation 
from WIMPs may be visible. A detection of annihilation products from multiple sources, 
in possible concert with detections from colliders and underground labs, will be required to 
conclusively establish the nature of the dark matter in the Universe. 

Each source of annihilation radiation has its advantages and disadvantages. Because of 
its close proximity and high dark matter density, the flux from annihilation radiation will be 
the largest in the direction of the Galactic center. However, uncertainties in the empirical 
determination of the central density profile and in the contamination from gamma- 

ray sources that are not of dark matter origin [5] may hinder the extraction of a dark matter 
signal from this region. Both of these systematics may be somewhat alleviated by searching 
for annihilation radiation at a few degrees offset from the Galactic center but even in this 
region a full analysis and understanding of the spectrum of the astrophysical backgrounds is 
required. 

The high mass-to-light ratios and the relative proximity of the dwarf spheroidal satellite 
galaxies (dSphs) of the Milky Way make them also excellent independent targ ets that have 
been widely considered in the literature fl E i a Q E B H UJ, \m- Their 
status as potential targets for indirect detection has become even more interesting recently, 
nven that the known number of satellites has more than doubled in the past few years 



171 118L 1191 1201 . 12111 . coupled with the discovery that all of these satellites share a common 



dark matter mass scale within their central 300 pc [|22l] . While the actual signal of gamma- 
ray flux from the Milky Way dSphs is smaller than the flux from the Galactic center, the 
astrophysical gamma-ray backgrounds tend to be reduced in the direction of these objects, 
as most of them are located at high Galactic latitudes. Additionally, these dSphs have low 
intrinsic emission from astrophysical gamma-ray sources; not only do they have mass-to-light 
ratios greater than ~ 100 in most cases, but all of them within ~ 400 kpc have strong upper 
limits on HI gas content [|23|] . 

The astrophysical contribution to the calculation of the annihilation flux from any source 
can generally be divided into two components: one component arising from the smooth halo 
contribution, which is proportional to the density squared distribution in the halo, and one 
component arising from the flux due to bound substructures within each of these halos. While 
it has long been recognized that the presence of substructure in dark matter halos can have a 
significant effect on the annihilation rate of dark matter particles theoretical calculations 
of this boost factor from substructure have varied by orders of magnitude because of large 
uncertainties in both the density profile of substructures and their distribution within the parent 
halos [0, 10, 25? ]. Numerical simulations are now reaching the necessary resolution to 



characterize the density distribution of substructures 11261 12711 . resolving substructure down 
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to four levels in the hierarchy [|27ll . These, the highest resolution simulations, model the 
evolution of the dark matter and do not include baryons. In principle, simulations of this 
kind can be used to provide accurate estimates for the expected boost factor only in objects 
with negligible baryonic mass fractions (e.g., dwarf galaxies). However, in regions like the 
Galactic center where baryons play an important role in the dynamics, baryonic effects, such 
as encounters of substructures with stars in the Galactic disk or bulg e [28.1 > or the backreaction 
of the dark matter distribution in response to disk formation [|29|], must be included. Future 
simulations that include more physics and even higher resolution will be required to fully 
characterize boost factors in general circumstances. 

It is clear that, in order to extract an unambiguous detection of dark matter annihilation 
radiation, a full understanding of all astrophysical uncertainties is required. With a goal of 
characterizing these uncertainties, in this paper we present an algorithm for the calculation of 
the gamma-ray flux from dwarf satellites, accounting for both uncertainties in the smooth halo 
distribution and the halo substructure distribution. We introduce a method for scanning the 
parameter space and determining the best fitting dark matter distributions from the kinematics 
of stars in these satellites using a Monte Carlo Markov Chain (MCMC) analysis. We combine 
the parameter constraints from the satellite stellar kinematics with the constraints on the 
parameters of the Constrained Minimal Supersymmetric Standard Model (CMSSM), with 
a goal of obtaining predictions accounting in a realistic way for all relevant sources of 
uncertainty of the flux from the dwarf satellites. The predicted regions that we delineate will 
provide guidance for future gamma-rays experiments for testing the predictions of neutralino 
dark matter in the CMSSM self-consistently within the context of CDM. As examples, we 
apply our algorithm to two particular satellites, both of which are known to be strongly dark 
matter dominated: the classical satellite Draco, which is located at 80 kpc from the Sun, and 
Segue 1, which is a newly-discovered satellite at 23 kpc from the Sun. We provide these 
two as examples and leave the work of ranking all the satellites in terms of their expected 
gamma-ray luminosity to a future paper. 

As a part of our analysis, we provide a new analytic solution to the equation governing 
the boost factor from halo substructure. Our solution is particularly useful because it allows 
for mass functions and halo concentrations to be free functions of host halo mass that can be 
implemented at each level in the substructure hierarchy. This is particularly important when 
implementing results from recent numerical simulations which show that the normalization 
of the mass function is reduced by up to 50% at the next level of hierarchy [|27ll . Using this 
solution, we show that the uncertainty in the boost factor is dominated by the extrapolation of 
the dark matter halo concentration versus mass relation down to mass scales that are currently 
unable to be resolved in CDM simulations. Assuming an optimistic power law extrapolation, 
we find mean boost factors ~ 20, in agreement with recent numerical extrapolations for Milky 
Way mass halos [27]. Assuming a concentration-mass relation that is linked to the small-scale 
power spectrum as in the model Bullock et al. [30], however, leads to boost factors of order 
unity. As an additional application of this analytic formula, we solve for the minimum dark 
matter halo mass at each point in the CMSSM parameter space, and find that the typical range 
for the minimum mass CDM halo is 10"^ - 10"^ Mq. This result updates and generalizes 
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the prev ious calculations of the minimum mass halo in the context of Supersymmetric CDM 

iQSHHIii. 

The paper is organized as follows. In Section[2lwe review the formalism for determining 
the gamma-ray flux from dark matter annihilations. In Section [3j we review our assumptions 
for the CMSSM and the SuperBayeS software that scans the CMSSM parameter space. 
In Section m we review the formalism to determine the best fitting dwarf satellite dark matter 
halo profiles, under the assumption of the CDM model. In Section[5]we present our calculation 
of the probability distribution for the boost factor and the resulting differences relative to the 
smooth halo flux. In Section |6] we discuss some detection prospects for present observatories, 
and in Section|7]we present our conclusions. Our main conclusions for the detectability of the 
flux are summarized in Figure [TT] 



2. Annihilation Flux and Gamma-ray spectrum 

Following standard methods IH [s^ the gamma-ray flux from particle annihilations can be 
derived via 

< cry > (E) r^'=^max n n(+ 

o(£) = — g^^/ da' J dQfR[e'-e) plMVmW{e), (i) 

where € is the line-of-sight distance, £+ = J cos 6 + sjr^^j^ - d^ sin^ 0, d is the distance to the 
galaxy, 6 is the line-of-sight angle from the center of the galaxy, < cry > is the average cross 
section for annihilation, rf OM is the tidal radius of the dark matter halo, and is the WIMP 
mass. Here, 

f"^ dN 

NJE)= —dE (2) 



is the number of photons above energy E produced per annihilation and the resolution window 
function is 



^/-\ In 2 



(3) 



re 

For Fermi, Ores is approximately 10 arcminutes at the energies we consider. It is clear 
from Equation ([T]) that an accurate prediction for the flux entails incorporating relevant 
uncertainties both for the astrophysical quantities and for the particle physics model. In the 
present analysis we constrain the density of the dark matter halo, p{l), from the kinematics of 
the stars, while the annihilation cross section, WIMP mass and annihilation channels we use 
are derived in the following section. 



3. The Constrained Minimal Supersymmetric Standard Model 

For the calculations in this paper, we will assume that the dark matter particle consists of 
the lightest stable supersymmetric particle. More specifically, we focus on the case of the 
CMSSM. In this section, we review the basic definitions and parameters of the CMSSM, 
discussing specifically the SuperBayeS code that we use to explore the CMSSM parameter 
space and how this feeds into our calculation of the flux. 
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3.1. Neutralino Dark Matter in the CMSSM 

Super symmetry (SUSY) provides a compelling and well-motivated extension to the Standard 
Model that naturally contains a dark matter candidate [36]. Supersymmetry postulates a 
symmetry between bosons and fermions - every boson (e.g. gauge bosons) has a fermonic 
partner (e.g. gauginos) and every fermion (leptons, quarks) has a bosonic partner (sleptons, 
squarks). The particularly well-studied R-parity conserving, weak-scale softly broken SUSY 
models provide both natural solutions to the "fine tuning problem" and a natural dark matter 



candidate [|37l] . The former is achieved through cancellation of quadratic divergences of 
one loop quantum corrections to the Higgs mass. The latter is due to a conserved discrete 
symmetry (R-parity) that prohibits the lightest SUSY particle from decaying to only SM 
particles. 

Enlarging the particle sector in this manner greatly increases the number of free 
parameters that specify the model; even the most minimal form of SUSY (MSSM) introduces 
over a hundred new parameters. Such a large number of free parameters makes the 
efficient exploration of the MSSM parameter space very challenging. The naive method 
of exploring the likelihood surface on a regularly-spaced grid is clearly inadequate, as the 
required computational effort scales exponentially with the number of dimensions considered. 
Furthermore, present-day constraints on SUSY phenomenology are fairly indirect and do not 
allow for meaningful constraints on models with so many degrees of freedom. A popular and 
well motivated simplification is achieved by demanding SUSY parameter unification at GUT 
(Grand Unified Theory) scales. This so-called constrained MSSM (CMSSM) now limits the 



number of SUSY parameters to 4 continuous and 1 discrete parameter 113811 : the common 
gaugino mass, mi/2, common mass for scalars, mo, trilinear scalar couplings, Aq, (all of which 
are specified at the GUT scale, Mqut - 2 x lO^^GeV) the ratio of the higgs expection values. 



tanyS, and the sign of the "/i term", sgn(//) 1371 1381]. We shall denote the CMSSM parameters 



as 

^ = {mo, mi/2, Ao, tan/?, sgn(ju)} . (4) 



It has been recently demonstrated [|39l . |40L l4ll . |42L 1431] that the value of some SM parameters 



can strongly affect the predictions for some of the observable quantities, in particular the relic 



neutralino abundance (see Fig. 4 in Roszkowski et al. [|41l]). Therefore, it is not sufficient to 
specify the values of ^ and fix the relevant SM parameters to their current best-fit values, 
but the latter must be introduced as "nuisance parameters" and marginalized over, in order to 
account for their impact in the predictions. The most relevant SM parameters are 

^ = {mj , mh, as , a^m) , (5) 

namely, the top quark mass, the bottom quark mass, the strong coupling constant and 
electromagnetic coupling constant, respectively. 

In the context of the CMSSM, specification of the parameter set C^, ^) allows 
for derivation of a full suite of predictions for low-energy observables. The 



package SuperBayeS , developed by Ruiz de Austri et al. [|40|] and Trotta et al. 11431] embeds 



several codes in a MCMC framework to derive from (^, ^) the SUSY mass spectrum (via 
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SoftSusy [44]), the neutralino relic abundance (via DarkSusy 1451] or MicrOMEG As ll46ll47ll ). 



50 



SUSY corrections to various Higgs sector quantities (employing FeynHiggs [|48l. 149 
and branching ratios of rare processes (using Bdecay [52]). The CMSSM and SM parameter 
space is explored by SuperBayeS using an MCMC Metropolis-Hastings algorithm, or, more 
recently, by employing the more efficient and robust "nested sampling" algorithm [|43Ll53L l54l]. 
The parameters are then constrained by applying all available constraints on the low-energy 
observables, including the WMAP 5-years determination of the relic abundance, sparticle 
and Higgs mass limits, branching ratios of rare processes, electroweak observables and direct 



constraints on the SM quantities (for a detailed discussion of the likelihood, see 11401 141114311) 



3.2. Priors in the CMSSM 

The final outcome of the CMSSM analysis is a list of samples drawn from the posterior 
distribution, P(^, ^\D) obtained via Bayes' theorem: 

/'(^, ^ ID) oc X(^, ^)P{^, ,jV) , (6) 

where D denotes the combined data described above, XC^, -yV) = P{D\^, ^) is the 
likelihood function and P{^, ^) is the prior distribution for the CMSSM and SM parameters. 
From the posterior one can then derive the probability distribution for any function of the 
quantities .yV) one is interested in, for example the neutralino-proton interaction cross- 
section (relevant for direct dark matter detection experiments, see 115 511 ). the gamma-ray and 
antimatter flux from the galactic center (relevant for indirect detection searches i56|] ) and 
Higgs-sector physics (of interest for Higgs-boson searches [|5^ ). 

The posterior in Equation ([6]) should be dominated by the likelihood, X, so that the prior 
influences vanish for strongly constraining data (for more details, see e.g. [i58il ). However, it 
has been found that this is not currently the case for the CMSSM — namely, the available data 
are not sufficiently constraining to determine the CMSSM parameters in a prior-independent 
way (see Trotta et al. [|43n for a detailed analysis). This means that some of the posterior 
constraints on ^ are somewhat dependent on the chosen prior distribution. The fundamental 
reason for this is that the mapping between high-energy CMSSM parameters and low energy 
observable quantities is highly non-linear, due to the nature of the Renormalization Group 
Equations, and therefore even fairly strong low-energy constraints are relatively mildly 
informative about the quantities one is interested in, namely ^. It is however expected that 
this issue will be resolved once the LHC will deliver direct observations of the SUSY mass 



spectrum ll59ll . 



3.3. CMSSM samples and derived quantities 

For the results in this paper, we use a nested sampling chain for the parameter space spanned 
by C^, ^), containing approximately 45, 000 samples. We assume throughout a positive 
sgn(yu) (motivated by consistency with the measured anomalous magnetic moment of the 
muon). We adopt the chain resulting from the analysis in Trotta et al. 14311 . with a flat prior 
with the limits denoted in table [D (the "4 TeV" limits in Ruiz de Austri et al. koll ). 
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From those CMSSM chains, we derive for each sample the value of the WIMP mass, m^, 
its annihilation cross section, < cry >, and number of photons produced in the annihilation 
above 1 GeV, Ny{l GeV). As discussed below, we choose 1 GeV because it provides a 
conservative lower bound for the expected signal energy window from CMSSM dark matter 
annihilation. Furthermore, we compute the value of the minimum mass halo, mmin, as 
explained in section [54l below. 



4. Dwarf Satellite Kinematics 



We follow the standard formalism for analyzing stellar line-of-sight velocities from dwarf 
satellites. In this section, we review the relevant formulae so as to establish notation and 
conventions that we use throughout this paper. 



4.1. Theoretical Modeling 



We consider the satellites as two-component systems consisting of stars and dark matter 
(e.g. [|60|, l6l|]). The potential is assumed to be spherically symmetric, and the system is 
taken to be completely pressure-supported (no rotational support). This is seen to be a good 
description of many of the dwarf satellites [1621 . |63L 16411 . With these assumptions the Jeans 
equation is 



dip-^crl) 



GM{r) 



(7) 



dr r 

where is the stellar density, cr^. is the stellar radial velocity dispersion, /3 = I - crj/crj is the 
velocity anisotropy parameter, and cr, is the stellar tangential velocity dispersion. The mass of 
the system, M(r), is defined as the total dynamical mass, which is the sum of the contributions 
from the stars and the dark matter. The line-of-sight velocity dispersion is 



o-th 



f 

Jr 



idr 



(8) 



where R is the projected radius onto the sky and /*(i?) is the stellar surface density. The use 
of subscript "f/z" will become apparent below when Equation ([8]) is fed into our statistical 
analysis. 

Observationally, the integrated masses within ~ 300 pc for all of the Milky Way satellites 
are very similar, consistent with ~ 10' Mq independent of the dwarf galaxy luminosity [1220. 
To first order, this fact simplifies the selection of the best targets for flux detection to those 
that are closest to the Sun [ |65ll . However, as we discuss below, including the effects of a 
prior relation for the maximum circular velocity distribution for dark matter halos somewhat 
modifies this simple estimate. In general, the primary factors that determine the best sources 
are those objects that have the best signal-to-noise ratio, accounting for the astrophysical 
backgrounds. Some often-discussed targets include Segue 1 (23 kpc), Ursa Major II (32 
kpc), Willman 1 (38 kpc), and Coma Berenices (44 kpc) (e.g. [|65l. l66l l67n) all of which were 



discovered since 2004. The half-light radii of these objects are ~ 10-100 pc, and given their 
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high velocity dispersions of ~ 4 - 6 km/s moA . these objects are consistent with being dark- 



matter dominated dSphs 11221. 167l l68h. The nearest of the more well-known (classical) dSphs 



include Ursa Minor (66 kpc) and Draco (80 kpc); previous calculations show that these two 
objects have similar predicted fluxes (l^. 

Since our main goal in this paper is to discuss the methodology for robustly predicting 
fluxes and including the boost in the calculation of the gamma-ray flux, we restrict our analysis 
to two example dSphs: Draco and Segue 1 . The kinematic constraints we derive for Draco are 
much stronger, as this object has been well-studied both from the standpoint of its photometry 
and kinematics. Segue 1 is a more-recently discovered satellite that appears to be strongly 
dominated by dark matter and the least luminous galaxy known [|670. However, there are only 
24 stars with measured velocities in Segue 1 [|67ll . and, as we show below, the errors on its 
mass and flux are much larger than the respective values for Draco. The surface densities of 
these objects are fit by King [69] and Plummer [70] profiles, respectively. For Draco, the King 



core radius is ryti„g = 0.18 kpc, and the King tidal radius is = 0.93 kpc [jyill . For Segue 1, the 
one-component Plummer fit gives a Plummer core radius of r^/ = 35 pc [|72|l . In the Plummer 



profile. falls ofl" as 1 /r^ in the outer regions, so there is no natural definition of the stellar 
tidal radius, in contrast with the King profile. In this case, we conservatively assume that the 
stellar tidal radius is given by the position of the outermost observed star, which is located at 
a projected radius of 7? = 50 pc 1671] . In principle, the stellar surface density parameters could 
also be estimated jointly with our other model parameters, however in the present work we 
choose to fix them instead to their best fit values as given above. 

To model the density profile of the respective dark matter halos of these objects, we use 
an Einasto profile, which is defined as 



p(r) = p, exp \ -In 



- -1 



(9) 



The Einasto profile has been shown to be a good fit to CDM halos with the Einasto index, 
n, ranging from ~ 3 - 7 [|73L l74l l75|l . For our purposes, this profile is also convenient for 
two separate reasons: (1) the profile has a well-defined mass, which is important when we 
calculate substructure boost factors below, and (2) the profile does not diverge towards the 
center of the halo, which is convenient when calculating the gamma-ray flux. As presented 
in Equation Q, there are three parameters in the Einasto profile that we must determine 
from the data: the log-slope index n, the scale radius r^, and the scale density, p^. It will be 
convenient in the following to replace the Einasto density and radius variables (p„ f^) with 
the implied halo maximum circular velocity and radius of maximum circular velocity (Vmax, 
Vax)- The dark matter halo density profile is then specified in terms of the parameter set 



^yif — {n, Vrnaxj 



(10) 



The final quantity we must specify is the velocity anisotropy, which enters both directly 
in the Jeans equation in Equation (|7]) and in the equation that relates the observed line-of- 
sight velocity dispersion to the underlying stellar radial velocity dispersion Equation ([8]). This 
quantity is unconstrained by line-of-sight velocity data uM, so in order to allow for general 
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models, we model the velocity anisotropy as 

with four free parameters, 

r = {/3o,fioo,rf,,r]}. 



(11) 



(12) 



We note that this parametrization is slightly more general than that used in Strigari et al. [|6ll] . 
allowing here for the power law index rj in addition to the anisotropy scale radius, rg, and the 
asymptotic inner and outer anisotropics, /3q and /3oo. It should also be noted that there is an 
intrinsic degeneracy between the logarithmic slope of the density profile in Equation (HI) and 
anisotropy. Even in the simplified case of constant anisotropy, this degeneracy restricts how 
well the central slope n of the halo may be contrained from stellar kinematics 1761]. 



4.2. Likelihood function ofdSph model parameters 

Above we have discussed our theoretical modeling; we now turn to our description of the 
line-of-sight velocity data [1671. 17711. We begin by noting that the line-of-sight velocities from 
dSphs are well-described by Gaussian distributions [62]. The observed velocity distribution 
is a convolution of the intrinsic velocity distribution, arising from the distribution function, 
and the measurement uncertainty from an individual star. The probability of obtaining a set 
of data d given a set of model parameters J^, Y is described by the likelihood [llOn 



1 



exp 



1 (di-uf 
Icrl . + (^2 . 

th,i m,i 



(13) 



where are the parameters that describe the density profile of dark matter and Y are the 
stellar velocity anisotropy parameters. The product is over the set of n stars, and u is the 
bulk velocity of the galaxy in the direction of the observer. As expected, the total error at a 
projected position is a sum in quadrature of the theoretical intrinsic dispersion, o-th,i(J^, y), 
and the measurement error cr^,,. 

Often, kinematic data from dSphs is presented in terms of the velocity dispersion in 
bins of projected radius. In these cases, it is useful to have an expression for the likelihood 
function similar to Equation (fT3l) that is free of any terms associated with the measured 
velocities of individual stars. An expression of this form can be found by replacing u and 
cr by their respective maximum likelihood values, m and &, where the latter quantities are 
obtained from a standard maximum likelihood procedure using Equation (031 ). Proceeding 
with this approximation, and also neglecting the measurement uncertainty in comparison with 
the intrinsic dispersion (which is a good approximation for the bright satellites), the likelihood 
function in Equation (fT3l) can be reduced to 



'2 



(14) 
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logio(Mass(30pc) [Mq]) logi(,(Mass(300pc) [M^]) 

Figure 1. Posterior probability distribution for the mass within 30 pc for Segue \{left panel) 
and the mass within 300 pc for Segue 1 {right panel). The four curves in each panel 
assume different Bayesian priors: a uniform prior in V^^^^ (black, solid) V^^^ (red, dashed), 
y~!i^ (green, dot-dashed), and In(ymax) (blue, dotted). The prior distributions are truncated 
at Vmax - 3km s"' as described in the text. Increasing negative powers of Vmax causes 
the posterior to be more "biased" toward lower mass solutions. As a result, the posterior 
corresponding to these different priors differ. 



This is now a product over the number of bins in projected radius, A'^, each with a velocity 
dispersion &b. The number of stars in a given bin is A^^,. (Note the difference in the 
normalization relative to the expression presented in Strigari et al. [65], due to a typographical 
error in Strigari et al. ^6^). Note that the quantity cr,k,i = cr,k,i{M', Y) is computed 
via Equation ([8]). As we did above in the context of the CMSSM parameters, we now go over 
to the posterior probability distribution function (pdf) for the parameters of interest, which is 
again given by B ayes' theorem 

where P{Jf), PiY) are the prior pdf's for the halo and velocity anisotropy parameters, 
respectively, which we take here to be uncorrelated. We will deal with the issue of priors 
in more detail below. 

The task is now to explore numerically the posterior. Equation (fT5l) . in order to determine 
credible regions on the model parameters. In previous similar work involving parameter 



estimation from dwarf satellites [|61L I65I]. the likelihood function was directly integrated over 
the range of model parameter space. This method was accurate but could be time-consuming, 
particularly in the case of large numbers of parameters to marginalize over. Rather than direct 
numerical integrations, in this paper we explore the posterior distributions using Markov 
Chain Monte Carlo techniques. Before discussing our MCMC methodology, we now turn 



to the discussion of the priors entering Equation (|T5l) . 
4.3. Priors for dSph parameters 



Choosing a prior in accordance with the physical situation and our degree of prior beliefs on 
it is an important aspect of Bayesian analysis, as variations in the priors can lead to sizable 
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7 7.5 
logi(,(Mass(300pc) [Mg]) 

Figure 2. Posterior probability distribution for the mass within 300 pc for Draco. The 
four curves in each panel assume different Bayesian priors: a uniform prior in V^^^ (black, 
solid) (red, dashed), V^l^ (green, dot-dashed), and ln(yn,ax) (blue, dotted). The prior 
distributions are truncated at Vm^x - 3km s"' as described in the text. Note that the trend in 
mass within 300 pc with prior reverses here compared to the case in Figure[T] This traces back 
in part to the fact that the mass is best constrained near twice the half-light radius [78]. For 
Draco, 300 pc is within the half-light radius and when the CDM rmax- V'max prior is imposed, 
lower Vmax halos are forced, on average, to be more concentrated at 300 pc. In Figure [1] 300 
pc is beyond the half-light radius of Segue 1, and priors that favor larger Vmax give larger 
extrapolated masses. 



differences in the posterior whenever the likelihood is not very strongly peaked. We account 
for this prior information in our notation by the inclusion of the appropriate infinitesimals. 
For example, a "uniform" prior probability in x will be denoted as P{x) = d{x) whereas a 
uniform prior probability in ln(x) is represented as P{x) = d[ln{x)] = d(x)/x. Regardless 
of these definitions, if data strongly constrains some parameter, or a given combination of 
parameters, then the prior information should not have much bearing on the result, as the 
posterior is dominated by the likelihood. However as we will see for Segue 1, with only 24 
radial velocities spanning only ~ 50 pc in radius, priors can have a significant effect. 

Guidance as to how to choose the priors for our dSph model parameters can be gleaned 
from cold dark matter simulations which give precise predictions of halo abundances given a 
halo shape and mass. We note however, that these simulations do not provide the probability 
of observing a galaxy, which itself depends not only on dark matter physics, but also on 
star formation and baryonic physics. But if we assume that primarily the halo mass, and not 
shape, affects stellar physics we may draw the halo profile prior probabilities directly from 
simulations, with some additional simple inputs to account for gas physics. 

As discussed above, we describe the halo in terms of the parameters n, v^x. ^max- We 
model the r^^x prior (conditional on the value of Vmax) as a log-normal distribution, which 
provides a good description of the Vax - Vmax relation measured in the Via Lactea ll79ll and 
Aquarius 1271] simulations over the entire range of Vmax- From the Aquarius simulations pvl] . 
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we estimate this relation, given a Vmax, to be 

pr IV ^ / [log(w) - l-351og(yn,ax) + 1-75]^ ] 

^(Vaxl Vmax) °^ exp I — \ d{\og w), (16) 

I log(rmax) ) 

where cr\og(r,,,.^^) = 0.22 is a conservative scatter in the Aquarius subhalos for the entire range 
of Vmax- There are no statistics published (as of this paper) for the parameter n in the Einasto 
profile, though from Springel et al. jivll it is reasonable to assume a uniform prior in 1 /n [i.e. 
P{n) oc d{l/n)] limited to the range 0.1 < l/n < 0.5. 

The choice of prior for VVax is a more delicate issue, for a couple of reasons. One issue 
relates to the probability that a given CDM halo has a particular value of Vmaxj and a second 
issue relates to how well the line-of-sight velocity data itself constrain particular values of 
Vmax- On the former point, Vmax is the primary parameter that relates to the astrophysical 
process of low-mass galaxy formation: small galaxies with more shallow potential wells 
are expected to have low star formation rates, so the actual Vmax prior is expected to be 
more shallow than that inferred from the CDM subhalo mass function, P(Vmax) oc V^ax^^^max 
(derived from the relation N(> Vmax) oc V^l^ of Springel et al. lilTIl ). On the latter point, 
from the perspective of the line-of-sight data, large Vmax values are not constrained by the 
data. Given the form of the Vax prior, large y^ax values correspond to large Vax values that 
may fall outside the extent of the stellar profile. Therefore, these values cannot be directly 
constrained by data and accordingly become dominated by the prior. 

Figure \T\ and Figure [2] exemplifies this behavior. The left panel of Figure [U shows 
the posterior constraints on the mass of Segue 1 within 30 pc with four different prior 
assumptions: /'(V^ax) oc diV'lJ, diV'lJ, diV^l^), and dilnV,-^^,^). The right panel shows 
the mass within 300 pc using the same respective priors. Prior choice has little effect on M(30 
pc) since this value is well constrained by line-of-sight velocity data. In contrast the M(300 
pc) posterior is dominated by prior choice, as the radius of 300 pc lies outside the measured 
stellar distribution. In contrast Draco, whose stellar extent extends beyond 300 pc, M(300 
pc) is well constrained. We note that the prior behavior in Figure \T\ and Figure [2] does not 
contradict the results presented in Strigari et al. where a uniform prior in ln[A/(300pc)] 
was taken for the entire satellite population, and no CDM-motivated priors were considered. 
Rather Figure[T]and Figure |2] pertain to the specific case of the prior assumed in Vmax- 

The above situation does, however, present a dilemma when considering priors in the 
quantity Vmax^ the actual V^ax prior for the observable galaxy is likely to be more shallow 
than the prior that comes from the predicted CDM V^ax prior for all substructure (uniform 
in y^ax)' but more shallow priors will give more statistical weight to parts of V^ax parameter 
space not well constrained by data. Because of these factors, the best that can be achieved 
with line-of-sight data is a lower limit to the model parameters. Thus, in this paper we use 
the CDM prior /'(Vmax) oc V'^ax'^^max with an imposed cut-off of 3 km s~'. This cut-off seems 
reasonable, as ~ 3 km s"' is expected to be a conservative lower bound to the Vmax values 
below which gas is able to condense into halos [80]. As we show, the imposed cut-off of 3 
km s ' does not affect the predicted flux from Draco, in which the posterior is dominated by 
the line-of-sight data. However, for Segue 1 (with only 24 stars), the issue is more subtle, as 
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logio(Mass(300pc) [Mg]) 

Figure 3. Posterior probability distribution for the mass within 300 pc for Segue 1, 
assuming an isotropic velocity distribution (the magenta dot-dot-dot-dashed line) and the 
model in Equation (fTTl i (black solid line). Both curves assume a uniform prior in V^^^. Both 
curves assume a uniform prior in V^-^ with an imposed cutoff below Vmax = 3 km s ' . 



the CDM prior becomes more dominating with decreasing values of Vmax- This is would be 
particularly true if the Vmax cut-off extended down to arbitrarily low values, though our choice 
of a (physically-motivated) cut-off somewhat reduces the effect of the prior, even for Segue 1 . 
In summary, our prior for the CDM halo parameters is given by 

PiM') = P(r^JV^,.)P{V^,,)P(l/r]) (17) 

with P(rmaxlV'max) as in Equation (fT6l) and /'(Vmax), Pi^h) according to the prescriptions given 
above. 

Finally, we have little physical basis to choose the prior for the anisotropy parameters 
'f . Additionally, these parameters are not well constrained by the data. Fortunately, these 
parameters, unlike the halo parameters, do not have a direct effect on the derived flux or mass 
(only indirectly through the Jeans equation). Thus, we expect that the choice of these priors 
will not have as large of an impact on the result as compared to the halo priors. Figure [3] shows 
the difference between the Segue 1 mass withing 300 pc using our prior assumption and the 
mass assuming an isotropic velocity distribution. Although our prior assumption biases the 
probability distribution to lower mass values, the effect is not as extreme as with the V^ax 
prior. 

The prior choices for all of the parameters (including particle physics quantities) 
considered in this paper are shown in Table [B 

4.4. Monte Carlo Markov Chain methodology 

Here we review the MCMC formalism necessary for our analysis; we refer to the papers 
referenced for more details. 

The goal of an MCMC algorithm is to generate a series of points in parameter space 
(called "a chain") with the property that their density is distributed according to the posterior 
pdf one wishes to explore. Then from the chain of "accepted" points, the marginal probability 
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Table 1. Summary of model parameters and the priors imposed on them. Unless otherwise 
stated, the prior pdf is flat within the given range. 



CMSSM Parameters, 


Priors Assumed 


Notes 


mo 


50 GeV < mo < 4 TeV 


CMSSM scalar mass 


mil 


50GeV<mi/2 <4TeV 


CMSSM gaugino mass 


Ao 


lAol < 7 TeV 


Common Trilinear Coupling 


tan/? 


2 < tan/? < 62 


Ratio of Higgs vev 




SM Parameters, ./f 


Priors Assumed 


Notes 


M, 


160 GeV <M, < 190 GeV 


Top quark mass 


nib 


4 GeV <mt<5 GeV 


Bottom quark mass 




127.5 < < 128.5 


EM coupling const. 




0.10 < ff, < 0.13 


Strong coupling const. 




DM Halo Parameters, ^ 


Priors Assumed 


Notes 


n 


0.1 < 1/n <0.5 


Einasto index, see Equation (|9]l 


''max 


see Equation (fTSIl 


Used to derive rs 


V 

^ max 


flat in y,74; V^^^ > 3km/s 


Used to derive ps 




Anisotropy Parameters, f 


Priors Assumed 


Notes 




< 77 < 3 


see Equation (fTTI) 




0.01 kpc < < 100 kpc 


Anisotropy scale length 




-2 < A) < 


Central anisotropy 




-3 < ySoo < 1 


Outer anisotropy 



distribution for each of the parameters is recovered by simply binning the points in the chain, 
and ignoring the uninteresting coordinates (two-dimensional distributions are obtained in a 
similar manner). In our case, we wish to explore the joint parameter space spanned by the 
particle physics model parameters, ^ and and by the dwarfs model parameters, and 
So we are dealing with a total of 15 parameters. Because the CMSSM likelihood and the 
stellar kinematics likelihood are independent, the joint posterior factorizes as (more details 
are given below) 

/'(^, ^, jif, y\D, d) oc p(^, ,j^)P{,3e, r )X(^, ^)X(^, Y). (i 8) 

A great advantage of the MCMC procedure lies in it efficiency, for the computational effort 
scales roughly proportionally with the dimensionality of the parameter space being explored, 
rather than exponentially. The true power of MCMC methods, which we specifically utilize 
in this paper, lies in the fact that in addition to obtaining the distributions for the model 
parameters, the probability distribution for any function of the model parameters is obtained 
by simply determining the function at each of the accepted points in the chain. In this way, 
we may easily determine the distribution of any parameter that is derived from our base 
set of parameters by post-processing the chain of accepted points. Examples of derived 
parameters that we are interested in for the purposes of this paper include the dark matter 
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mass of the dwarfs or the gamma-ray flux. While the former quantity is just a function of 
the parameters that describe the dark matter halo, Y) the latter quantity requires us to 
combine the probability distributions on ^) with the probability distribution from the 
dwarf kinematics. To understand how these probabilities can be combined, and quantities 
such as the flux can be robustly calculated, we appeal to the general properties of conditional 
distribution functions. 

A Markov chain of the joint posterior distribution P{'^, jV , , d). Equation (fTSl) . 
can be obtained from the combination of the two posteriors Pi!^ , ^\D) and P(J^, y\d) as 
long as the two joint probability distribution can be factorized as in Equation (fTSl) above. This 
is, in fact, the case here since the likelihood for the particle physics model is unaff"ected by 
the stellar kinematic data. It is also the case that the halo density profile and stellar velocity 
dispersion anisotropy parameters are not affected by the particle physics constraints, hence 
the two likelihoods are independent. On the other hand, if we had included, for example, 
gamma-ray flux upper limits from Imaging Air Cerenkov Telescopes (ACTs) [81] in the 
particle physics likelihood, this would couple the two separate parameter spaces and invalidate 
the above decomposition. At present including such data is not necessary because the upper 
limits are well above the CMSSM parameter space, however it will be desirable and indeed 
necessary in the future. 

To obtain a Markov chain on the dSph model parameters, we opted for a combination of 
the slice sampling [ | 82|] and the Metropolis-Hastings [83] algorithm. The advantage of slice 
sampling is that, unlike the Metropolis-Hastings algorithm, its efficiency is not strongly linked 
to the proposal pdf. Whereas, with a good proposal pdf, the Metropolis-Hastings algorithm 
converges faster to the desired distribution. Thus, we obtain an initial proposal pdf using 
slice sampling and then employed it in the Metropolis-Hastings algorithm to derive our final 
pos terior. (For the actual slice sampling methodology, see Neal [|82ll and Lewis and Bridle 
[84]. For the Metropolis-Hastings methodology see Christensen et al. [|85ll . Lewis and Bridle 
ISfl], Baltz et al. IsTI], and Lewis and Bridle iMD for details). 

In practice, we found that our slice sampling run took 3-4 likelihood evaluations per 
point and offered fairly good convergence. However, our subsequent Metropolis-Hastings run 
had an acceptance rate of 30%-50% with excellent convergence. Nine chains of 30,000 points 
were obtained, thinned, and then combined with the SuperBayeS chains (as per the method 
outlined in Appendix A. 2 1. We refer to Appendix A.l for a more detailed discussion of the 
convergence criteria we use for our chains. 



5. Boost from Halo Substructure 



In the previous sections, we have outlined our modeling of the dwarf halos from stellar 
kinematics and our method for scanning the CMSSM parameter space and used this modeling 
to predict the flux under the assumption of a smooth halo. The final ingredient we must add 
to the flux predictions is the boost from halo substructure. The goal of this section is to derive 
the probability distribution for the boost factor, accounting in the most reasonable possible 
manner for the astrophysical and particle dark matter uncertainties that enter the calculation. 



Indirect Dark Matter Detection from Dwarf Satellites 



17 



5.1. Defining Boost 

Dark matter halos form hierarchically, and this results in a population of surviving 
gravitationally-bound substructure. High resolution dissipationless numerical simulations 
reveal substructures in z = Milky Way-size halos with a mass spectrum that rises towards 
smaller masses with A^(> m) oc nT" and a ~ 0.9, down to the smallest masses that can 



be resolved, currently m > 10~^Mhost ~ lO^M© [|27LlZ9|]. However, as we discuss below, 
this mass is some ten orders of magnitude larger than the minimum mass expected for 
CDM structure, mmin- While numerical simulations that focus on nested regions within very 
high redshift halos have demonstrated that halos with masses close to the filtering mass do 



survive the initial process of halo formation [|88L 1890 , more modeling will be required to 
better understand the survival probability, density structure, and precise mass spectrum of the 
smallest CDM substructures at z = 0. 

Because substructures themselves were assembled from smaller units prior to infall into 
their host, we expect a hierarchy of substructures within substructures that extends down to 
mmin- The first explorations of this hierarchy of mass functions and substructure distributions 
is just now becoming viable in state-of-the art simulations [27]. The dark matter halos of 
dSphs are substructures of halos of the Milky Way, therefore their substructure fractions 
and mass functions are less well explored, but there are clear expectations. Relative to 
substructures within the Milky Way halo, substructures within dSph halos are expected to be 
older, and this gives them more time to assimilate into the smooth component of their hosts. 
Moreover, the halos of dwarf satellites do not get replenished by accreted field halos to the 
extent that an isolated field halo would. Both of these effects will act to reduce substructure 
fractions in dSph halos relative to isolated galaxy halos. 

Following previous treatments, we define the boost B such that total gamma-ray flux O 
from a halo of mass M is 0(M) = [1 -I- B{M, m^in)] 6(M). Here 6 is the flux that comes from 
a smooth halo of mass M and and the boost includes a contribution from all subhalos larger 
than mmin, which is set by particle physics. We begin with the formulation in Strigari et al. 



jlOll . but note that the formalism of Fieri et al. [|250 provides similar results. Adopting the 



above definition, we may write the boost as an integral that accounts for substructures going 



down the CDM mass hierarchy [|lOl] : 



0(M)5(M, mn,in) = AM" [1 + B{m, m^^^)] (^im)m-'-"dm. (19) 

Here we have assumed that halos of mass M host substructures of mass distribution 
dN/dlnm = A{M/m)" for m < qM and (for now) we assume a self-similar substructure 
hierarchy. Written in this way, we can see that the total boost in gamma-ray luminosity 
depends sensitively on a competition between the smooth flux <i)(m), which tends to decrease 
towards smaller masses, and dN/dlnm, which rises to small masses. For halos described by 



NFW density profiles [e.g.l30|,LZ4|] with a concentration-mass relation that follows c(m) oc m 



with yu ~ 0.1, we expect 0(m) oc m^ with ^ ^ I - 2.2// [|lOl] . Setting ^ = 1, we find 
that the boost 5(M, mmin) is approximately proportional to (M/mmmT~^ ■ Since the ratio 
(M/mmin) ~ 10'^ for host halos of relevance, the precise value of the quantity a - ^ (~ 2.2//) 
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at the smallest mass scales becomes crucial in determining whether the boost is significant 
or negligible. By making the most optimistic assumptions possible {q = \, self-similar 
substructure hierarchies, and optimistic assumptions for the density structures of small halos 
with ju ~ 0.1) Strigari et al. [10] and Kuhlen et al. have shown that the boost should be 
no larger than about 100, and more typically < 10 for mmin ~ IO'^Mq. 



5.2. Analytic Solution for Boost 



The above solutions to Equation (fT9l ) [llOl 12611 were obtained by numerical integration. Here, 
we present an analytic solution to Equation (fT9l) . and discuss its utility. 

We begin by assuming that the flux at a fixed host halo mass scales roughly as a power 
law, Q){M) oc . With this, we may rewrite Equation (fT9l ) as 

DQ))' = A [q^-" exp{b) + D(ln q + b)] (20) 

where b = lnq + \nM - Inm^in and D(M) = M^-"B(M). Note that B{m^Jq) = D(0) = 0. It 
should be noted that B{mmm/q) = and not B{mmm) = is the appropriate boundary condition 
(as was assumed in Strigari et al. lEoil ). Using these boundary conditions Equation (|20l) can 
be recursively solved as: 

' M < m^Jq 



B(M,m^n) 



a4 



I qM \° 



1 



a 



1 

lq<M < m^Jq 



(21) 



This forms a set of functions, B{M,m^^^ = {5o(M, mmin), 5i(M, m,-nin), . . . , 5„(M, mmin)} 
where each function B^iM, m^in) is only valid in the interval m^i^lq'^ < M < m^i^lq'^'^^ . 

Conceptually, the 5,,'s represent the amount of substructure included in the calculation 
of the boost. For example, Bq{M, m^^^ represents the boost of a halo with no substructure, 
and thus we set Bq(M, m^in) = 0. Bi{M, nimm) represents the boost with the inclusion of only 
the subhalos whereas BiiM, m^nin) includes only the subhalos and sub-subhalos. B^iM, m^in) 
now can be related through Equation (|2T]) by 

D„iby = A [q^-" expib) + D„-i(ln q + b)]. (22) 

We may now solve for B^iM, m^in) by taking the Laplace transform of Equation (l22l) (see 



AppendixTB) ). After inversion, this yields (for n > 0) 

Aq^-" ' 



" 1 / 
BniM,mmin) = 2] (7TT)! ( 



a 



y z, i^-a) In 



q'M 



m„ 



(23) 



(24) 



where y{a, x) is the lower incomplete gamma function defined by 

y{a,x)= I x"'^ Qxp{-x)dx. 
Jo 

In the above analysis, we have assumed that the mass function is self-similar for all 
levels of substructure. However, this is unlikely to be the case and simulations [1271 . 17911 see 
less substructure in subhalos than in host halos. To account for the fact that the mass function 
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may differ for each level of substructure, it would be useful to perform the same analysis 
for the case where the mass function varies independently for each level of substructure. We 
define the mass function to be 

a mm \m I 

at level i. Here, / = would apply to the parent halo whereas / = 1 would apply to the 
subhalos and so on. Also, we let q qt following the same notation. Using the same analysis 
as above, we can rewrite Equation (1231) as 

,,(M,™.., = |]^(M^),(,tf-<,)l„(M)) (26) 



where A, = YVj=i ^iid qi = Y\'j=i 1j- For completeness, we present the a, = ^ solution : 

. i,.(ln(M.)y 

B„(M, m^in) = 2] (27) 



The utility of the analytic solution derived above stems from the fact that it lets us account 
for each level of substructure separately. We find that in the most extreme circumstance of 
^-a = -0.2 and m^nm = 10"'°Mq, inclusion of only the subhalos and sub-subhalos leads to an 
accuracy of 98.4%. Thus, it is not necessary to go below n = 2, for accurate boost predictions 
- for most part it is suflicient to just resolve subhalos to estimate the boost. 



5.3. Annihilation luminosity - mass relation 



In order to compute the boost, we need to estimate the gamma-ray luminosity from a subhalo 
of the halo under consideration. The luminosity depends strongly on the log-slope of the inner 
density profile of the subhalo and the eff"ect of the tidal forces of the parent halo. In order to 
estimate these effects, we consider a simplified model of subhalo distribution and tidal effects. 

We consider two possibilities for the spatial distribution of the subhalos. One is motivated 
by the current high-resolution simulations of Milky Way sized halos [1751 l90l l9lL 19211 . where 
it was found that the number density of subhalos with masses greater than ~ 10"^ the mass 
of the parent halo falls off as 1 /r^ at large radii. At small radii, tidal forces of the parent halo 
will destroy the subhalo and hence reduce the population in the inner most regions. However, 
it is not clear that the smallest subhalos extending down to earth mass subhalos should have 
this, more diffuse, spatial distribution. We therefore also consider the other extreme where 
the spatial distribution of the subhalos follows the smooth distribution. We choose the second 
possibility for the final computations because the boost is primary affected by the smallest 
mass halos. We randomly generate subhalos according to the chosen spatial distributions. 

Before it fell into the host halo, we allow for the log-slope of the subhalo inner density 
profile to range between a = 0.8 and a = 1.4, motivated [|93|l . We note, however, that more 
recent simulations indicate that there is no asymptotic inner slope [1271 l74l 19411 . The density 
profile assumed for the subhalo when it was outside the host halo is 



Pair) = po(r/ror"(l + r/ro) 



(28) 
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where a = I corresponds to the NFW profile and the concentration is defined by the ratio 
of the virial radius and the scale radius as c = R^iJrQ. We assume a power-law virial mass 
function for the subhalo before it fell into the host halo and generate masses randomly for the 
subhalos. Since we start by picking the virial mass of the halos, we use the observed mass 
function for field halos. However, there is a huge extrapolation involved - we assume that 
the power-law extends all the way down to the minimum halo mass where the mass function 
gets truncated. Given the virial mass, we find the concentration of the halo using the field 
concentration-mass relation including the scatter. While these relations were derived for the 
NFW profile, we adapt them for the profile in Equation (1281) . Assigning a concentration is 
the most uncertain part of the calculation; in fact, as we will see, the flux depends sensitively 
on the concentration of the smallest halos. 

Following the above procedure, we are able to set ro and po, and from this we then find 
the corresponding ergodic distribution function, fa(E). Given the implied tidal radius, we 
adopt the following simple model for tidal effects wherein any given tidal radius defines a 
relative energy Eq below which the distribution function of the subhalo drops to zero. This is 
analogous to the King models for the stars. That is, 

f(E) = f,(E)-iE>Eo, (29) 
= OV£<£o- 

(30) 

Two points are worth noting: (1) Eq is fully specified by the ratio of the tidal radius to ro, and 
(2) the product por^ has to be unchanged because the tidal effects do not change the innermost 
regions of the subhalo. We find that the density profiles resulting from the distribution function 
in Equation (l29l) are well fit by 

Psub(r)=p.(r)exp(-(r/r,)"), (31) 

where is defined here as the tidal radius and « is a function of Eq (or equivalently r,/ro). For 
any given po and ro, this model then provides a relation between r, and the tidally truncated 
mass of the subhalo. If r, is smaller than ro, we consider the subhalo destroyed. 

To fix the mass and tidal radius for a given po, ro and distance from the center of the host 
halo, we iteratively solve the Jacobi relation between rf and mass, under the approximation 
that both the host halo and the subhalo are isothermal. We allow for the fact that the satellite 
could have been closer in its orbit by choosing the closest distance to the center along the orbit 
to be a factor less than one, chosen randomly from to 1, times the current distance. 

This is now a fully specified model and we may now predict the gamma-ray annihilation 
luminosity as a function of the tidally truncated mass of the subhalo. At the high mass 
end of the subhalo mass function, the predictions will have large scatter because of sample 
variance. At smaller masses, where most of the flux comes from, the large number of subhalos 
contributing to the flux means that the predictions essentially have no scatter. 

As a means of illustrating the important uncertainties associated with subhalo density 
structures we explore two simple models for the concentration-mass relation. The first is 
a power law (PL) model c(m) oc m that is normalized to match the results of N-body 
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simulations over the mass ranges that have been probed by direct simulations (m > IO^Mq). 
The second is the simple analytic model Bullock et al. ll30ll with the specific implementation 
of Maccio' et al. [[9511 . This model (denoted BOl below) also matches simulations down to 
the mass-scales explored by simulations, but links the value of c to the spectrum of CDM 
density fluctuations via estimated collapse times E|. While the power-law extrapolation leads 
to concentrations of order 500 for earth-mass halos, the BOl model concentrations are much 
lower, and plateau for small masses because the slope of the power spectrum of density 
fluctuations is more shallow on small scales. 

We note that the flux depends primarily on the profile around r^ and smaller. The 
concentration sets po since po oc c^ and the flux is proportional to pl/^y We find that the 
scaling formula suggested by Strigari et al. [llOn . 0(M) oc Mc(M)^--^, works well (here M 
is the subhalo tidally-truncated mass) for power-law c(M) functions. We did not find any 
systematic differences in this scaling for the two assumptions about the spatial distribution 
of the subhalos. However, we do find that the more concentrated spatial distribution (NFW) 
implies systematically larger fluxes (by a factor of about 2). In the self-similar calculation 
below we only include the scaling information. 

We do not consider the effect of the assumptions about the spatial distribution of subhalos 
on the signal-to-noise. In particular, if the distribution of the subhalos is as shallow as 
l/(r^ + r)^, then much of the signal (even for more moderate boosts of order 1) will come 
from the "outer" regions. As we integrate the signal outwards from the stellar core, two effects 
must be considered. One, the signal-to-noise depends sensitively on the angular acceptance 
region around the satellite. If the background is mainly extragalactic (or at least constant 
across the angular region of the galaxy), it will scale linearly with the angular region covered. 
The scaling of the signal with the angular acceptance can only be estimated if the spatial 
distribution of the subhalos of a satellite is known. Two, we cannot accurately predict the 
tidal radius beyond which both the dark matter density profile of the satellite galaxy is cut-off. 
For the boost calculation in the MCMC exploration, we assume that the spatial distribution 
of the subhalos follows that of the smooth halo. We reiterate that the relevant quantity is the 
distribution of small subhalos that cannot yet be resolved by simulations and those subhalos 
are the ones that contribute dominantly to the gamma-ray flux. 



5.4. Minimum Mass CDM Halo 

We now have an expression for the boost as a function of (1) the host halo mass, (2) the 
mass function of CDM subhalos, (3) the concentration-mass relationship for CDM subhalos, 
(4) the minimum mass CDM halo, m^nm- The uncertainty in the fourth item arises from the 
unknown cut-off scale in the mass function of CDM substructure. As mentioned above, this 
cut-off scale is well below the resolution of modem day CDM numerical simulations of Milky 
Way like galaxies at z = 0. The smallest halo size is set by the horizon at kinetic decoupling 
in the early universe. The kinetic decoupling in turn depends on the scattering interactions 

t The choice of BOl is representative of many proposed analytic models that link c to the power spectrum. 
Among these it has the steepest faint-end slope. 
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of dark matter with standard model fermions, as well as the free- streaming after decoupling 

m^in = 5.72 X 10-2q^/j2^3/4 r^Vgj// j (32) 
Here geff = 10.75 is the effective degrees of freedom when the CDM particle freezes out, and 

2 

m 

C=^\M(t^O)\\ (33) 

Pan 

where \M(t 0)\ is the matrix element of elastic scattering in the limit of no momentum 
transfer (t ^ 0). 

Bertschinger Isil] has calculated mmin assuming a Bino-like neutralino, finding typical 
values mmin ~ 10""^ Mq for typical WIMP masses. We follow the kinetic decoupling 
calculation of Bertschinger [ISO but generalize it to include Wino and Higgsino-type 



neutralinos by determining C at each point in the SuperBayeS chain. In Appendix C 



we present the relevant Feynman rules and matrix elements that enter into calculating the 
scattering matrix in Equation (1331) . 

It is important to note that in our analysis of the scattering matrix element, we assume 
that before kinetic decoupling occurs, the WIMP interacts solely with electrons and neutrinos. 
In doing so, we have neglected the effects of muon scattering, which will be relevant if 
kinetic decoupling occurs at temperatures comparable to the muon mass. Including the non- 
negligible muon abundance would in turn modify Equation (|32l) . which is beyond the scope 
of our present analysis. 

In Figure m we show the resulting posterior distribution of Wmin accounting for the entire 
presently viable parameter space of the CMSSM. We find ~ 95% c.l. values for the minimum 
halo mass within the range ~ 10^^-10"^ M©. The features in the likelihood for nimm arise from 
the probability distributions in the CMSSM parameters; due to the non-linear transformation 
between these parameters and the minimum mass halo in Equation (|32|) any small features in 
the CMSSM parameters are strongly exacerbated in the m^in likelihood. 



5.5. Boost Predictions: Two models 

Running Monte Carlo simulations (see section [53l) . we found that the luminosity 0(M) has 
a power law behavior for both the PL and BOl models. In terms of the first-order estimates 
of the boost discussed in association with Equation (fT9l) above, the relevant combination of 
the mass-function log-slope (-a) and the luminosity function log-slope (f) takes the value 
a - ^ - 0.2 and 0.1, for the PL and BOl models, respectively. Figure [5] shows the boost for a 
10^ Mq halo across its relevant range for both models. The ensuing posterior pdf for the boost 
for Segue 1 for these two models is shown in Figure [6l The difference in boost between these 
two models is an order of magnitude and this underscores the large effect that the (uncertain) 
concentration of the lowest mass halos has on the overall boost. This leads to a natural 
uncertainty for any flux prediction and emphasizes the need to measure the concentration 
and mass function well for halos of the smallest size. In order to be conservative, for the 
remainder of this paper we assume or - ^ = 0. 1 . 
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Figure 4. Posterior probability for m,nin in the CMSSM, assuming flat priors and employing 
all available constraints. 
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Figure 5. Left: Halo concentration versus halo mass for our PL model (solid black line) and 
the BOl model (dashed red line). Right: The boost for an M = lO^M© halo that results from the 
PL model (solid black) (using Equation (|23] |) and the BOl (dashed red) concentration models. 
Note that both of these concentration-mass models are consistent with current simulations. 



6. Flux Predictions and Detection Prospects 

In the analysis above we have assembled the necessary ingredients to make gamma-ray flux 
predictions given constraints from the halo kinematics and the CMSSM. The corresponding 
kinematic constraints on the (smooth) dark matter halos of Draco and Segue 1 are illustrated 
in FigurelTl Now, as a final step in the process of making the flux predictions, we must specify 
the angular region around the center of the dSph within which the flux is calculated, .i.e. ^max 
in Equation ([T]). Determining the most optimal value of ^^ax which maximizes the signal and 
minimizes the background contribution depends on the experiment and backgrounds under 



consideration. Satellite experiments such as Fermi [|97d and ACTs, such as HESS 1981 ] and 



MAGIC [|99ll. are expected to have an angular resolution of ~ 0.2° and 0.1°, respectively. If 



the extent of the gamma-ray emission region for the dSphs is larger than these point- spread 
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Figure 6. Left: Posterior distribution value for the boost factor, assuming the concentration- 
mass model of BOl. Right: Posterior distribution value for the boost factor, assuming the PL 
concentration-mass model. 
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Figure 7. Allowed parameter space for the scale radius and scale density for the Einasto profile 
in Equation (|9]l for Segue I {Left) and Draco {Right). Inner and outer contours represent the 
68% and 95% c.l. regions respectively. 

functions the dSphs would be resolved as point sources. 

Our approach in this paper is to provide an algorithm for determining the optimal value 
of ^maxj given the best fitting halo parameters and the measured background spectra. More 
specifically, we determine the optimal value of ^^ax from our Markov chains, along with 
the estimation of the detection significance, <t = NJ y/N^ + Nb, where A^;, is the number of 
background photons and A^^ is the number of signal photons within a given angular acceptance. 
We note that the definition of detection significance here is only meant to be an approximation; 
more realistically the maximum signal-to-noise will depend on the angular distribution of 
the signal and background, and detailed detector specifications. Nonetheless, proceeding 
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^max [degrees] [degrees] 

Figure 8. The posterior probability distribution for the angular region within which the signal- 
to-noise ratio is maximized. Left panel is for Segue 1, and Right panel is for Draco. The sharp 
cut-off at low 0max is a result of the assumed point-spread function. 



with the above definition, at each point in our chains, we calculate the value of 6'niax for 
which the significance cr is maximized, and then construct the pdf of the quantity, ^max- 
In order to convert the flux to a detected number of photons, we take the total exposure, 
defined as the orbit-averaged efi'ective area times the observation time, to be 3 x lO^'cm^s ^ 
2000cm2 X 5years H- 

For our input background spectrum we perform a standard analysis utilizing the EGRET 
diffuse background measurements at high Galactic latitude. We consider just the diffuse 
extragalactic background, though note that including contributions from the residual Galactic 
component may increase the total background by flux about an order of magnitude. The 
diffuse backgrounds over the relevant energy range of ~ 1 - 100 GeV for WIMP annihilation 
will of course be mapped with even greater precision in the very near future with Fermi. 
The d iffus e extragalactic background is seen to fall off according to the power law dN/dE ~ 
E~^-^ [lod]; we integrate this background spectrum for energies > 1 GeV, and compare to the 
number of photons produced by dark matter over the same energy range. The resulting pdf 
for 0max is shown in Figure [81 As Figure [8] shows, accounting for the diffuse backgrounds, the 
optimal angular extent for each galaxy is ~ 0.2° (the PSF of Fermi), with a gradual tail that 
extends to larger angles. 

Given the distribution of S^ax for each galaxy, in Figure |9] and Figure \T0\ we present our 
results for the flux of gamma-rays from each galaxy with energies greater than 1 GeV. In 
these figures, we show the gamma-ray flux versus two CMSSM parameters, the neutralino 
mass and the spin-independent cross section. For these figures we assume the "BOl" boost 
model, where the minimum value for the halo mass is self-consistently calculated from the 
CMSSM parameter space. Because it has more line-of-sight velocities, the flux distribution 
for Draco is more strongly-constrained. As is seen, the predicted fluxes are similar for both 
galaxies, despite the fact that Segue 1 is more than a factor three times closer than Draco. 
The main reason that these fluxes are similar traces back to the prior used in Equation (fT6l) : 
given the velocity dispersion and half-light radius. Segue 1 prefers to reside in a halo with 
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Figure 9. The E > \ GeV gamma-ray flux versus the neutrahno mass for Segue 1 {Left) and 
Draco (Right). Both figures include the conservative boost model corresponding to the Bullock 
et al. 13011 halo concentration model (discussed in Section 15.5b and CDM prior (discussed 
in Section l431 l. Inner and outer contours represent the 68% and 95% c.l. regions respectively. 
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Figure 10. The E > \ GeV gamma-ray flux versus the spin-independent cross section 
for Segue 1 {Left) and Draco {Right). Both figures include the conservative boost model 
corresponding to the Bullock et al. |30] halo concentration model (discussed in Section l53b 
and CDM prior (discussed in Section |43] |. Inner and outer contours represent the 68% and 
95% c.l. regions respectively. 



lower Vmax relative to Draco. The corresponding one-dimensional flux distributions are show 
in Figure nn 

Using the definition of significance, a (defined above), and an assumed exposure, we 
can obtain a rough estimate of detection prospects. As an example, we consider the detection 
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Figure 11. Probability density for the flux of £ > 1 GeV gamma-rays from Segue I {Left) 
and Draco (Right), after marginalizing over all parameters of the CMSSM and the dark matter 
halos. The solid and dashed curves correspond to the boost model with the Bullock et al. 
1 301] and power-law halo concentration models (discussed in Section 15.51 1 respectively. The 
dashed-dotted line denotes the flux above which the signal to noise is about 3, using the 
estimated definition of signal-to-noise is section |6l Notice the double peaked feature in the 
Draco posterior. This is due to the fact that the astrophysical component of the flux is very well 
constrained by data causing the uncertainty in the flux to be dominated by the non-Gaussian 
CMSSM parameter space. The uncertainty in the flux for Segue 1, though, is dominated by 
the astrophysics and accordingly is very Gaussian in shape. 



prospects at a signal-to-noise greater than 3 (cr = 3) in five years (same exposure as given 
previously). Comparing to our flux predictions in Figure [HI we see that with the standard 
BOl boost model, the parameter space spanned by the flux posterior of Draco lies below this 
estimated flux limit (vertical dot-dash line) and only a small fraction of the parameter space 
reaches this limit for Segue 1. On the other hand, taking the optimistic PL boost model, in 
which the mean boost is ~ 20, both Segue 1 and Draco exhibit significant regions of their 
parameter space that are detectable given our estimated Fermi sensitivity. Specifically, for the 
PL boost model, we find ~ 20% of the Draco flux parameter space is detectable with signal- 
to-noise greater than 3 with five years of Fermi observation (and similarly ~ 13% of the Segue 
1 parameter space at the same level). We note again that the assumed priors (see below) have a 
large eff"ect on the predicted flux and our general strategy in this work has been to quantify the 
minimum expected flux. We also note that for our PL model we assumed a spatial distribution 
for the subhalos that tracks the underlying host halo mass distribution. Making this spatial 
distribution more extended will lower the flux from the central region. Our results for the BOl 
concentration boost model are broadly consistent with the results of Pieri et al. [fl6l], who also 
considered annihilation flux from Draco but with an energy threshold of 100 MeV. 

For the PL boost model, it would be interesting to also consider the flux of gamma- 
rays from annihilations in the unresolved substructure of the Milky Way and ensure that this 
is not in conflict with the measured EGRET background [|25l] . However, such an exercise 
is hampered by the facts that the Galactic disk could significantly change the substructure 
distribution as well as abundance, and there is no systematic way of taking this into account 
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at the present time. 

Of course, to extract a dark matter annihilation signal above background, a detailed 
comparison between all of the input spectra is required, rather than simply counting photons 
above some energy threshold. In order to make a detailed comparison between spectra, it will 
be important to determine the characterstic energy for photons from dark matter annihilation. 
We define the characteristic energy as E^^^, the energy at which the quantity E^dN/dE peaks. 
For each point in our Markov chain we determine the characteristic energy, and the resulting 
pdf for Ejaax is given in Figure [T2l 

Finally, we re-vist our above discussion of priors in the context of flux predictions. We 
see from Figure [H] that when the astrophysical data is well constraining the majority of the 
flux uncertainty comes from the particle physics parameter space. This severely reduces the 
impact of the astrophysical priors on the overall uncertainty of the flux. But for galaxies 
like Segue 1 with only 24 line-of-sight velocities, the likelihood for its parameters is much 
less constraining and thus astrophysical priors have a large impact on the flux posterior. 
In Figure [T3j we show how the assumed prior aff"ects the resulting calculation of the flux 
from Segue 1 and Draco. The curves in Figure \T3\ are for the same priors as in Figure \T\ 
It is thus clear that, given the small number of stars from Segue 1, the assumed prior has a 
significant effect on the flux calculation. The CDM V^ax prior forces Segue 1 to reside in 
halos with smaller V^ax and this results in smaller predicted fluxes. We note that the physics 
of faint-end galaxy formation should change this prior significantly and would generally push 
the prior towards favoring larger Vmax values. Thus our chosen priors that have the effect of 
providing the minimum expected flux from Segue 1 . Future data sets for this object, and all 
other ultra-faint satellites, will be crucial for constraining the dark matter mass and estimating 
the gamma-ray fluxes. 

The priors on the CMSSM parameter space may also play an important role in 
determining the flux expectations from the dwarfs [I43I]. We have not undertaken a systematic 
study of the effect of these priors on the fluxes here. We chose a uniform prior in all the 
CMSSM parameters including the masses. 
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Figure 13. Effect of the prior on the predicted flux for Segue 1 (Left) and Draco (Right). 
Again, a uniform prior in V^l^ (black, soHd) (red, dashed), V^l^ (green, dot-dashed), and 
In(Vinax) (blue, dotted) is assumed. As in the case of the mass (see Figure [T]i, flat priors in 
increasing negative powers of y^ax causes the Segue 1 posterior to be more biased than the 
kinematically well-constrained case of Draco. For these figures, we have set the boost factor 
equal to zero. 



7. Discussion and Conclusions 

Constraining the properties of the dark matter particle in indirect detection experiments 
will require a firm understanding of the astrophysical uncertainties that contribute to the 
flux. Dwarf satellites of the Milky Way are particularly interesting targets in this regards 
because they are the most dark matter dominated objects known and they are largely free 
from astrophysical uncertainties which result from the presence of baryonic physics. In this 
paper, we have taken an step towards quantifying uncertainties in flux predictions from dSphs 
by providing a framework within which both particle physics and astrophysics uncertainties 
can be included at once. 

We have combined our MCMC method for determining the dark matter distributions 
of dwarf satellites using stellar kinematics with the SuperBayeS MCMC package which 
determines the preferred ranges for parameters of the CMSSM. We have focused on two 
specific dSphs, Segue 1 and Draco, as example cases. Our methods allow us to provide a 
broad outline for the prospects of detection of these satellites with gamma-ray experiments 
such as Fermi. 

The main results of our paper can be summarized as follows: 

• We find that both Draco (at 80 kpc) and Segue 1 (at 23 kpc) are expected to have grossly 
similar fluxes, though the flux from Segue 1 is subject to larger uncertainties because of 
its relative lack of kinematic data and smaller stellar extent. For the most conservative 
assumptions. Segue 1 prefers to reside in a halo with lower V^ax relative to Draco, and 
therefore it has a similar overall flux despite its relative proximity. However, for a flat 
prior in logCVmax), the flux from Segue 1 can be much larger than Draco (see Figure [T3]). 
This result motivates future observations of stellar kinematic data for Segue 1. We note 
that our results for the flux, unless otherwise stated, are based on the more conservative 
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prior and conservative boost model, and hence they quantify the minimum expected flux 
for the assumed CMSSM priors. 

• We have provided the first self-consistent calculation for the boost in the flux signal from 
halo substructure, taking into account both the CMSSM model and the recent results 
from high resolution numerical simulations. We show that the dominant uncertainty in 
the boost calculation comes from the assumed halo concentration versus mass relation 
for halo substructure on mass scales down to the scale of the minimum mass halo. If 
we assume a model that links halo concentrations to the power spectrum, then we obtain 
typical boost factors of order unity. If we instead assume a power law continuation of the 
concentration-mass relation down to the minimum mass halo, the average boost factor 
increases by an order of magnitude to ~ 20 (see Figure (6]). This boost would be reduced 
if the spatial distribution of the smallest subhalos is more diffuse that the smooth dark 
matter halo component, which is what we have assumed in this paper. We also note that 
our analytic solution for the boost shows that resolving subhalos, and in some cases sub- 
subhalos, is sufficient to get an accurate estimate of the boost. These facts motivate future 
high resolution simulations of halo substructure that can measure the concentrations of 
small dark matter halos directly, as well as map their spatial distribution within the host. 

• We have provided a broad outline of the prospects for detection of these satellites with 
gamma-ray experiments, focusing specifically on Fermi. We find that, given the diffuse 
backgrounds, the most optimal solid angle to view these galaxies is ~ 0.2 - 0.3° (see 
Figure [H). Optimistic fluxes for these galaxies are approximately a few times 10"^' 
cm"^ s~'. For the boost model with the power-law extrapolation of the concentration- 
mass relation and a subhalo spatial distribution that tracks the underlying host mass 
distribution, we estimate a ~ 20% chance for a ~ 3(T dark matter gamma ray signal 
from Draco after 5 years of observation. This expectation is highlighted in Figure El 
We note that given the observed uniformity in the central density of the dark matter halos 
of the dwarfs [|22|]. it should be possible to stack them and increase the signal-to-noise. 

• We have provided an updated calculation for the minimum halo mass in CMSSM. Our 
results broadly agree with those of Profumo et al. [|34l] . and are typically lower than the 
masses in Bertschinger i35|] who only considered bino-like neutralinos. We find that the 
minimum mass CDM halo lies in the range 10"^ - 10"^ Mq (see Figure (6]). The inclusion 
of this effect results in slightly larger boost estimates. 

The methods presented here provide a concrete methodology for addressing the 
uncertainties inherent in dark matter indirect detection from both particle physics and 
astrophysics perspective. The focus on Segue 1 and Draco was meant to illustrate how two 
vastly different kinematic data affect the predictions for the gamma-ray flux from dark matter 
annihilations. It is very interesting to note that Fermi observations of the dwarf satellites of 
the Milky Way could be relevant for constraining the Super symmetric parameter space. 
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Figure Al. Power spectra normalized to a variance of one of a Gaussian random deviate (left), 
and of our converged chain for the parameter rp (right). Notice for small k (large scales) that 
the MCMC chain becomes flat like that of the true random deviate. The dashed lines are the 
best fit of Equation dA.lb . 
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Appendix A. MCMC methods and tests 

Appendix A. 1. Convergence Test 

In this appendix we briefly discuss our criteria for convergence of our Markov Chains. 

Points close to each other in a Monte Carlo Markov Chain are highly correlated and thus 
do not represent the distribution being explored. But, as the chain becomes larger, points far 
from each other become less correlated and eventually become random deviates of the target 
distribution. When this happens, the ch ain is said to have "converged." This eff"ect can be 



seen in the power spectrum of the chain [|l01|l . The power spectrum of a truly uncorrected 
random sample will be flat whereas a random walk power sectrum has the form P{k) oc 
where k is inverse of the distance between the points. Thus, for a converged MCMC chain, 
the power spectum will be constant at small k but will vary as k^ at large k (see Figure IaT)) . 
For a converged chain, the power spectrum can be fitted by the function 

{k^lkf 

P(k) = Po — . (A.l) 

Using the best fit values from equation [A. li the conditions for convergence for a chain of size 
were set to be 

(i) /^fc*^>20 (A.2) 

2n 



Indirect Dark Matter Detection from Dwarf Satellites 



32 



(ii) ^ ^ ^ ^-^^ ^^-^^ 

where the power spectrum is normalized to a var iance of one and Equation (lA.ll) is fitted over 
the range <k < lOk* (see Dunkley et al. lIlOlll for details). It should be noted however that 
this test does not guarantee that the chain has sampled the entire parameter space, but only that 
its autocorrelation lengths is sufficiently short to ensure independence among the samples. 



Appendix A.2. Combination Methodology 

If two parameter sets are completely independent, then the statistics of the combined 
parameter set is independent of whether each parameter set is derived independently or 
jointly. Thus, if a posterior is separable (!P(0, ifj) = P,/,(^)P^(>ff)), then the total posterior 
may be found by the combination of the individual distributions. But care must be taken 
combining chains with different weights. If the resulting combined weights differ from 
the individual distribution weights, then the resulting combined distribution runs the risk of 
differing from the original, pre-combined, distributions when marginalized over. A solution 
to this problem is to explicitly repeat points in each chain (thereby separating the weights 
out among the repeated points) in such a way that the two chains sets of weights coinside. 
For example, lets consider the two chains = {1,2,3,4} and i/f = {a,b,c,d} with weights 
"IV^ = {4,3,1,1} and "IV^ = {1,2,3,3} respectively. These chains are equivalent to the 
lengthier chains of (p = {1,1,1,2,2,3,4} and i/^ = {a, b, c, c, d, d, d} with the common weight 
set of "IV^ ,/, = {1, 2, 1, 2, 1, 1, 1}. Now these two chains can be combined into the merged 
chain of {(1, a), (1, Z?), (1, c), (2, c), (2, J), (3, J), (4, J)} with the above mentioned weight set. 
This method will only, at maximum, double the chain length while preserving the individual 
parameter distributions. 



Appendix B. Boost Derivation 



In this appendix, we derive Equation (|23]) . We start by taking the Laplace transform 
of Equation (|22)) : 



sFn(s) 



+ Aq'Fn-i(s) 



(B.l) 



where F^is) is the Laplace transform of D^ib). By recursively substituting F,j{s), we can 
derive 



Fn(s) = 
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(B.2) 



There are two poles in Equation (|B.2I) - one at = and one at s = ^-a. The inverse Laplace 
transform is the residue at both these poles. After a bit of algebra, this becomes 
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(B.3) 
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Figure CI. Tree level feynman diagrams for the s, u, and t channels respectively for the 
process x + ^ ~* A" + ^■ 



where d\ represents the i-th derivative with respect to s. The first and second terms in Equation 
(IB.3I) come from the the residues at * = and s = ^ - a respectively. Equation (IB.3I) in its 
current form is very impractical considering for each term i- \ derivatives must be performed. 
This motivates finding a simpler relation for the first term in Equation (IB. 31) . Letting c = ^ - or 
andy = q'M/m^in, 
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This leads to 
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7h,(^-a)ln' 



(B.4) 



(B.5) 



The argument in the incomplete gamma function can either be positive or negative and follows 
the recursion relations 



y{i,c\ny) = (z - l)y(/ - 1, clny) - (clnj)' ^y' 
y{\,c\ny) = 1 -y''. 



(B.6) 
(B.7) 



Appendix C. Scattering Matrix Calculation 



In this appendix we derive the cross section for neutralino-lepton elastic scattering. We 
closely follow Chen et al. fl02]'s procedure in finding the neutralino-neutrino scattering 
matrix element with the exception that we generalized their result for any lepton. We start 
with the amplitudes fo r s, u, and t channels repr esent ed in Figure ICll (see Jungman et al. 
Bal . Haber and Kane llOj], Gates and Kowalski u04\ for the appropriate Feynman rules). 
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The t channels contributions due to an Higgs boson exchange become negligible with the 
assumption of highly relativistic leptons (e.g. the coupling constants vanish for 0). The 
matrix amplitudes due to each individual diagram are 

= 'J^^' + <^"3 "[c-i (x;. + w;. U2 (c. 1) 
Mu= - i—^u, (x;. + w^;. .p,) u,u, (x',,.,Pl + w;,^-„p«) «2 (c.2) 

w in - 



Mr = i 



2 cos 6„ t 



m% 



where 



Pl = 
Pr = 



1 -y5 
2 

1+7^ 



Cl = 



Cr = 







-u^r (clPl + crPr) U2U^y, (O^j^Pl + 0';„^Pr) u, (C.3) 

(C.4) 
(C.5) 
(C.6) 



(C.7) 



for Ve,V^, Vr 

1 for e, IX, T 

for Ve,V^, Vr 

2 sin^ 6„ for e, fi, t. 

and the coupling constants X'^ij„, ^'ujn^ ^nd are defined in reference L36]. Here, the 

indices t, i, j, and n represent the family type (charged lepton or neutrino), the lepton flavor, 
the slepton flavor, and the scattering neutralino. Here, C = y^y" is the charge conjugation 
matrix. Using the above amplitudes, the combined matrix element is 



where the individual components to the matrix element are 
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In the limit of highly relativistic leptons and by assuming » E[, the specific momentum 
scalar products can be approximated by 

2 



{pi ■ p^) - s - u - 

(Pi ■ Pi) - iP2 ■ p-i) - iP3 ■ Pa) - (Pi ■ Pa) - m^Ee 

t 2 
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This leads to a matrix element of 
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